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ABSTRACT 

Euler's equations of motion are derived exactly for a rigid, triaxial, internally 
frictionless neutron star spinning down electromagnetically in vacuo. It is shown that 
the star precesses, but not freely: its regular precession relative to the principal axes of 
inertia couples to the component of the radiation torque associated with the near-zone 
radiation fields and is modified into an anharmonic wobble. The wobble period t\ 
typically satisfies t\ <^ 10 _2 to, where tq is the braking time-scale; the wobble amplitude 
evolves towards a constant non-zero value, oscillates, or decreases to zero, depending 
on the degree of oblateness or prolateness of the star and its initial spin state; and the 
(negative) angular frequency derivative Cj oscillates as well, exhibiting quasi-periodic 
spikes for triaxial stars of a particular figure. In light of these properties, a young, Crab- 
like pulsar ought to display fractional changes of order unity in the space of a few years 
in its pulse profile, magnetic inclination angle, and Co. Such changes are not observed, 
implying that the wobble is damped rapidly by internal friction, if its amplitude is 
initially large upon crystallization of the stellar crust. If the friction is localized in 
the inner and outer crusts, the thermal luminosity of the neutron star increases by a 
minimum amount AL ~ 3 x 10 31 (e/10 _12 )(a;/10 3 rads _1 ) 2 (T<}/l yr) _1 ergs -1 , where 
e is the ellipticity and is the damping time-scale, with the actual value of AL 
determined in part by the thermal conduction time r con d. The increased luminosity is 
potentially detectable as thermal X-rays lasting for a time ~ max(Td, r con( j) following 
crystallization of the crust. 
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1. INTRODUCTION 

Does the angular momentum vector of an isolated neutron star change orientation as the 
star spins down? Early attempts to answer this question focused on the evolution of the angle a 
between the star's rotation and magnetic axes, which can be measured from pulsar polarization 
swings. Davis &: Goldstein (1970) showed that a tends towards zero on the braking time-scale 
if the star is a rigid sphere or a fluid body in hydrostatic equilibrium — an unrealistic scenario 
which leaves all except the youngest pulsars as aligned rotators, contrary to observation. Goldreich 
(1970) observed that the crystalline crust of a neutron star supports shear stresses, thereby 
preventing a fraction of the hydrostatic bulge from aligning with the instantaneous rotation axis 
and establishing non-hydrostatic differences between the principal moments of inertia. In the 
absence of internal friction, such a triaxial star precesses about its principal axis with a period 
that is short compared to the braking time-scale, and the (fixed) angle between the principal 
axis and magnetic axis determines whether the precession amplitude (related to a) increases or 
decreases under the action of the braking torque (Goldreich 1970). 

In reality, a neutron star is not internally frictionless. Elastic strain energy is dissipated in 
the crust as the non-hydrostatic deformation discussed above migrates around the star while it 
precesses (Goldreich 1970; Chau & Henriksen 1971; Macy 1974). There is also dissipation due to 
imperfect coupling between the differentially rotating crust and superfluid core (Shaham 1977; 
Alpar & Ogelman 1987; Link, Epstein Sz Baym 1993; Sedrakian, Wasserman & Cordes 1998). Both 
types of friction damp any precession that is initially present on a time-scale that is, theoretically, 
much shorter than the precession period — a plausible explanation for why the pulse profiles and 
polarization characteristics of young pulsars do not change secularly over several years as one 
expects if the precession is undamped. The only isolated neutron star unambiguously known to 
precess is PSR B1913+16 (Weisberg, Romani k, Taylor 1989), where general relativistic effects 
are responsible. Tentative reports also exist of oscillatory variations in Crab and Vela timing 
residuals (e.g. Lyne, Pritchard & Smith 1988; McCulloch et al. 1990; Cadez, Galicic & Calvani 
1997), and some authors have inferred changes in a on the braking time-scale from braking-index 
measurements (Allen &; Horvath 1997; Link & Epstein 1997) and the evolution of pulsar radio 
beam statistics (Tauris & Manchester 1998). 

The radiation torque acting on the neutron star has been treated in an incomplete fashion in 
all of the above work. The torque has two components: the familiar braking torque, responsible 
for the secular spin-down of the star, and a component associated with the inertia of the near-zone 
radiation fields, sometimes misleadingly termed the 'anomalous torque' (Good & Ng 1985), whose 
effect is to make a spherical star precess about its magnetic axis. In all analyses to date, either 
the total radiation torque has been set to zero, in order to study free precession (Macy 1974; 
Shaham 1977; Alpar & Ogelman 1987), or else the near-field component has been neglected, in the 
belief that it exerts no significant influence on the rotation (Davis & Goldstein 1970; Goldreich 
1970; Chau & Henriksen 1971; Sedrakian et al. 1998). A careful study of both torque components 
in the context of a spherical star was carried out by Good & Ng (1985) for a magnetic dipole, 
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a magnetic quadrupole, and a hypothetical distribution of magnetospheric currents. Casini & 
Montemayor (1998) recently explored some effects of the near-field torque on a composite body 
with a spherically symmetric crust coupled to a spherical core (cf. de Campli 1980). 

In this paper, we demonstrate that the near-field component of the radiation torque strongly 
influences the rotation of an internally frictionless neutron star. In Section |^, we derive and solve 
Euler's equations of motion for a rigid, triaxial magnet and show that the regular precession 
relative to the principal axes of inertia couples to the near-field torque, causing the star to wobble 
anharmonically. Certain potentially observable properties of the wobble are explored in Section 
H, including its period and amplitude and the slow evolution of a and the angular frequency 
derivative to. The results are applied to pulsar timing and polarization observations, and to the 
internal structure of young neutron stars, in Section [||. 



2. ROTATION OF A RIGID, TRIAXIAL NEUTRON STAR 



In this section, the rotation of a rigid, triaxial body with an embedded magnetic dipole 



is treated analytically. Euler's equations of motion are written down in Section 2.1. Two key 
elements of the motion — the torque-driven radiative precession, and the inertial free precession 



are discussed in Sections |2.2| and 2.3 respectively and their time-scales are identified. Euler's 



equations are then solved approximately in Section 2.4 by a time-averaging technique for the 
special case of a biaxial star. The similarities and differences between this treatment and previous 



work are noted in Section 2.5. 



2.1. Euler's Equations 

Consider a rigid, triaxial star with principal axes ei, &2 and e3, corresponding principal 
moments of inertia I±, I2 and I3, ellipticities e = (I3 — and e' = (I2 — h)/ h, and average 

radius tq. We assume that the stellar magnetic field is dipolar and fixed in the star, and we restrict 
the magnetic axis m to lie in the plane spanned by ei and e3, at an angle x to 63; this entails a 
slight loss of generality. We also assume that the star is internally frictionless and that it rotates in 
vacuo, so that it is unaffected by magnetospheric currents; refer to Goldreich (1970) and Melatos 
(1997) for justifications of the latter assumption. It is important to keep in mind that, although 
we idealize the star as a rigid body for simplicity, in reality the ellipsoid of inertia is determined 



by an equilibrium between elastic and hydromagnetic forces, as discussed in Section 2.3 



In Appendix [A], we evaluate the radiation torque acting on the star from the electromagnetic 
fields generated by a magnetized, conducting sphere rotating in vacuo (Deutsch 1955; Melatos 
1997); it is assumed that the star's triaxiality can be neglected when calculating the fields and 
torque. Upon resolving the torque into components along the principal axes, we arrive at Euler's 
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equations of motion, 

(l + e')« 2 
(l + e)«3 

with 



(e' - e)-u 2 u 3 + (w r ) ^osx u 2 .F(:co)(-ni cos x + u 3 sinx) 
+ uG(x )u 2 (ui sinx + U3COSX) , 



-l 



-U 2 F(X )U2 



+ uG(xo) (-1*1 cos x + U3 sin x) («i sin x + u 3 cos x) 
-e'u\U2 - (woTo) -1 sinx u 2 F(x )(-ui cosx + ^3 sinx) 
+ uG(x )u2 (ui sin x + u 3 cos x) 



i^o) 
G(x ) = 



Xn 



b(x b - 3xg + 36) 
3(x§ + 6) 



+ 



3(xg + 1 
3 



2xg 



(1) 
(2) 
(3) 

(4) 
(5) 



3^ + 36) 15xoK + l) 

In equations (||) to @, u;o is the magnitude of the angular velocity vector w(t) at time t = 0, we 
define dimensionless variables u = u/lvq and xo = (roCt>o/c)ii, an overdot denotes differentiation 
with respect to the dimensionless time coordinate s = uot, and 



is the characteristic braking time-scale at t 
magnetic poles (Deutsch 1955). 



(6) 



0, in terms of the magnetic field strength Bq at the 



The form factors F{x$) and G(xq) reflect the structure of the near-zone radiation fields 
through their dependences on xq. The familiar braking torque, which is responsible for the secular 
spin-down of the star, is associated with terms proportional to F(xq), whereas the near-field 
torque discussed in Section [l] is associated with terms proportional to G{xq). In most applications, 
ro is taken to be the stellar radius R (cf. Section |3^ and Kaburaki 1981, Melatos 1997), yielding 
xq <C 1, F(xo) = 1/3 and G(xq) = 3/lCteo- The near-field torque is therefore much greater than 
the braking torque in this regime and acts on a commensurately shorter time-scale. Equations (Q) 
and (||) differ from the expressions F(xq) = 1/3 and G(xq) = 1/2xq found in previous works (Davis 
h Goldstein 1970; Goldreich 1970), partly because the treatment in this paper is not restricted to 
xq <C 1, and partly because we model the star's internal magnetization in a slightly different way, 
as explained in Appendix [A]. 



2.2. Radiative Precession 



The near-field component of the radiation torque causes the star to precess and nutate about 
its magnetic axis. We call this motion 'radiative precession'. To understand its origin, consider the 
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simple special case of a spherical star (e = e' = 0). In the regime xq <C 1, where the form factor 
uG(xo) fa 3c/10roaJo is independent of s, Euler's equations ([[]) to (||) have the exact solution (cf. 
Davis & Goldstein 1970) 



U 2 
U 3 



|exp 


~2i4 fi F(x )s~ 
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r 1/2 

> cos 
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(7) 

(8) 
(9) 



for the initial conditions u\ = Ui t o, ui = 0, 113 = u^fi 7^ at s = 0.0 A spherical star therefore 
precesses harmonically about m with period 



n 



( 207rr u;o \ 
\ 3n 3)0 c J 



TO < T , 



(10) 



and the precession amplitude decays exponentially on the time-scale ^tq/u^ q, as w aligns with m. 
Taking ro = R = 10 km (cf. Section 3.6) and I\ = 1 x 10 38 kgm 2 , we obtain numerically 



n = 4 x 10 13 

which is to be compared with 



( B 



r = 5 x 10 
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(11) 



(12) 



Radiative precession occurs due to the asymmetric inertia of the near-zone radiation fields of 
a rotating magnetic dipole (Goldreich 1970): the electromagnetic energy density £ is greater at 
the magnetic poles {£ fa Bq/2^q) than at the magnetic equator {£ fa Bq/8/Uq), translating into a 
fractional distortion e md fa (5£ / c 2 )^/ 1\ of the moment of inertia about m. The effect depends 
only on the radiation fields outside the star, which in turn are determined completely by the 
magnetic field at the stellar surface and the property that the star is a good conductor (Deutsch 
1955). The magnetic field inside the star does not influence the radiative precession. 

The near-field torque contributes terms proportional to U2U3 and U1U3 in 

d) and (|) 

respectively and therefore adds to, or subtracts from, similar terms arising from material 
distortions (i.e. e and e'). Thus a biaxial star (e' = 0) with x = and e = e rad does not precess at 
all, because the terms (e' — e)u2Us and eu\U3 are cancelled out exactly by the near-field torque. 



4 See Davis & Goldstein (1970) for a discussion of the singular special case 113,0 = 0. 
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2.3. Free Precession 

In the absence of the radiation torque, a triaxial neutron star precesses relative to its principal 
axes of inertia on a time-scale T2 = 2Tr/eu, where e is the non- hydrostatic ellipticity. A variety of 
non-hydrostatic mechanisms, many with geological analogues (Lambeck 1980), combine to deform 
the stellar mass distribution. We concentrate on elastic and magnetic deformations in this paper. 

The crystalline stellar crust supports shear stresses which prevent a fraction of the hydrostatic 
bulge from aligning with the instantaneous rotation axis. For a crust with uniform shear modulus 
[i, one finds e cr = 5flio 2 R 3 /4(l + p)GM, with p, = 38vr/ii2 4 /3GM 2 , where M is the mass of the 
star (Goldreich 1970; Lambeck 1980, p. 42). This elastic deformation yields a precession period 

t " = 1x1 ° 14 (io4^)"(t^Ft)" s < 13 > 

for neutron star parameters. If the crust is structured as a Coulomb lattice, one has jj, <^ 10 -4 and 
the non-hydrostatic fraction of the bulge is small; however, the exact value of jl is uncertain and 
may be much less than this upper bound. 

The magnetic field inside the star creates an additional deformation because non-radial 
field gradients (e.g. between the poles and equator if the field is a dipole) support non-radial 
matter-density gradients in hydromagnetic equilibrium. (This is not related in any way to the 
electromagnetic inertia of the external radiation fields discussed in Section |2.2j . ) The geometry 
and magnitude of the deformation is difficult to estimate, because little is known about the 
structure of the internal magnetic field. Thompson & Duncan (1993) argued that, if the internal 
field is generated after collapse in a convective dynamo, then it must be organized into randomly 
oriented loops ~ 1 km in size, each with field strength ~ 10 11 T; the neutron star rotates too slowly 
(Rossby number Ro ^ 10) to establish a coherent toroidal field at the base of the convection 
zone (cf. the Solar dynamo). On the other hand, if the internal field is generated before collapse 
in the progenitor star (e.g. in the convective outer envelope or hydrogen-burning core), a large 
toroidal field can grow (Ro <^ 0.1). Blandford, Applegate & Hernquist (1983) examined toroidal 
field generation by a thermoelectric dynamo. Several authors have raised the possibility of a very 
strong internal field (;> 10 10 T 3> Bq) in diverse contexts, including off-centred-dipole theories of 
the pulsar death line (Arons 1998), thermally regulated resurrection of a buried field (Muslimov & 
Page 1996), Ohmic decay in an anisotropically conducting core (Haensel, Urpin & Yakovlev 1990), 
crust-core coupling in Vela glitches (Abney, Epstein & Olinto 1996; cf. Easson 1979), and the effect 
of a virial field (10 14 T) on modified Urea cooling (Yuan & Zhang 1998) and the quark-hadron 
equation of state (Pal, Bandyopadhyay & Chakrabarty 1998). 

Assuming that the internal magnetic field is at least as strong as the surface field and roughly 
dipolar, we find that the hydromagnetic deformation satisfies e mag ;> e rad c 2 /c 2 , where c s is the 
isothermal sound speed (= 3 _1 / 2 c in a relativistic star); in other words, e mag exceeds e rad at least 
by a factor of order the ratio of the stellar to Schwarzschild radii. This yields an upper bound on 
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the precession period r™ ag associated with the hydromagnetic deformation given by 

r 2 maS < 1 x 10 13 ( 2 ( ~* s . (14) 

Note that the combined elasto-hydromagnetic deformation is triaxial in general. Biaxiality (e' = 0) 
is a good approximation only when one has e cr <C e mag (or else e cr ^> e mag ) and the internal 
magnetic field (or crust) is symmetric about a unique axis. 

Equations (|TT|), (|D|) and (O) reveal that, as a rule, the periods of the radiative precession 
and free precession are comparable. One has n ~ , if the internal and surface magnetic 
fields are of similar magnitude, and sometimes t\ ~ t 2 as well, e.g. for a one-second pulsar with 
/i = 10 28 Nm~ 2 and Bq = 4 x 10 8 T. In general, therefore, the 'free' precession is not free at 
all; rather, it couples to the near-field component of the radiation torque, and its character is 
modified significantly as a result. This coupling has been overlooked in the literature on aspherical 
rotators to date (Goldreich 1970; Chau & Henriksen 1971; Macy 1974; Shaham 1977; Alpar & 
Ogelman 1987; Sedrakian et al. 1998), and the remainder of this paper is devoted to exploring 
its consequences. It can only be neglected under certain circumstances, e.g. when the internal 
magnetic field is strong and the free precession is therefore fast (r™ ag ^ 0.05ri; see Section |3.2| ). 



2.4. Separating the Braking and Precession Time-Scales: An Approximate Solution 

of Euler's Equations 

Euler's equations (|l|) to (|3|) can be solved approximately by averaging over the precession 
period, exploiting the fact that both t\ and t 2 are small compared to tq (except in the regime 
xo ~ 1; see Section [3.6p , The analysis parallels that by Goldreich (1970), with one crucial 
difference: we account fully for the radiative precession in what follows, whereas Goldreich (1970) 
artificially suppressed it by averaging the braking torque and near-field torque over the free 
precession period. 

Let us restrict attention to a biaxial star (e' = 0) for the sake of simplicity. When the 'slow' 
braking terms in Euler's equations, proportional to (u>oTq)~ 1 F(xo), are neglected relative to 'fast' 
terms, proportional to e and (u!qto)~ 1 uG(xo), the equations reduce to the zeroth-order system 



Ul 

with 



e/3 1 u 2 [ui sin x cos x - (P - cos 2 x) u 3], (15) 
e P~ l \{ u l -u 2 )sinxcosx + (/? - cos 2 x + sin 2 x)u\Uz], (16) 
-eP~ l u 2 {ui sinx + u 3 cosx) sinx, (17) 



<i8> 



The parameter (3 is independent of s in the regime xq <C 1 and is positive or negative according 
to whether the star is oblate or prolate. Upon multiplying (15), (|l~6|) and (|i~7|) by ui, u 2 and U3 
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respectively, we find that 

TJ = u\ + u\ + u\ (19) 

is a constant of the motion; in other words, the total angular momentum of the star is unaffected 
by the near-field component of the radiation torque. Furthermore, dividing ( |l5| ) by (17) yields a 
second constant of the motion 



7 = (uisinx + u 3 cosx)^ 



, 



(20) 



which loosely measures the difference between the precession and magnetic inclination angles. 
One can solve (|l6|), (|i~9| ) and (|20|) simultaneously to obtain a first-order differential equation for 
U2, solvable by quadrature, whose solution is exactly periodic but anharmonic in general. For the 
illustrative special case x = 90° > /? > calculated in Appendix H, we find 



U-2 



U3 



—z cn 

v-j 



sn 



1/2 



dn 



(3 

y-j 



1/2 



S + $ 



1/2 



S + $ 



1/2 



S + $ 



(21) 
(22) 
(23) 



where sn, cn and dn are Jacobian elliptic functions with modulus k 2 = (7 + t](3)/(rj — 7)/?. Note 
that the phase $ of the oscillation is a third constant of the motion, related to u±, 112 and U3 in a 
complicated way. 

When terms proportional to (cjqTq) F(xo) are restored to Euler's equations, the above 
solutions remain approximately valid, but the constants of the motion are converted into slowly 
varying functions of s, viz., rj = Tj(s), 7 = j(s) and $ = <J?(s). The slow variation occurs on the 
braking time-scale tq, which is long compared to 77 ~ T2 ~ e^ 1 . Analytic expressions for 77, 7 and 
$ are obtained by substituting (|2"T|), ( p2| ) and ( p3| ) into (|l]), (^) and (|3|) and averaging over 
Results for x = 90° are given for reference in Appendix [b|. 



2.5. Comparison with Previous Work 

Equations to (|3|) reduce to the equations solved by Davis & Goldstein (1970) for the 
special case of a spherical star (e = e' = 0, and hence x = without loss of generality), and to 
those solved by Goldreich (1970) in the regime T2 <C t\ where the free precession is much faster 
than the radiative precession. The latter regime corresponds to a large elasto-hydromagnetic 
deformation, such as when the internal magnetic field is much stronger than the surface field (see 



Section 2^3), but it is not fully general. In contrast, the solutions presented in Section 2A and the 
numerical results in Section 3 address the general problem where T1/T2 is arbitrary, including the 
regime 77 ~ T2 where the free and radiative precessions couple together. Previous authors have 
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discussed the physical origins of t\ and T2 (Goldreich 1970; Chau & Henriksen 1971; de Campli 
1980). 

Our analysis cannot be compared directly with previous work treating the star as anelastic 
(e.g. Macy 1974) or with models featuring a core and crust coupled together (Shaham 1977; de 
Campli 1980; Alpar & Ogelman 1987; Casini & Montemayor 1998; Sedrakian et al. 1998). 



3. PROPERTIES OF THE ROTATION 

The motion analyzed in Section Q is characterized by several properties of potential 
observational significance which we now investigate, including the precession period (Section |3.1| ) 
and amplitude (Section |3.2j ), the evolution of the magnetic inclination (Section p.3| ), the different 
behaviour of oblate and prolate stars (Section |3.4| ), the phenomenon of 'pseudo-glitches' in Co 
(Section |3.5[) , and the effect of a corotating magnetosphere (Section |3,6| ). In what follows, we 
assume that the orientation of u is arbitrary at the time when the neutron star first crystallizes 
into an object with a rigid crust spinning down electromagnetically. In other words, 111,0, ti2,o 
and 1*3^ are assumed to be comparable; we do not choose a privileged initial spin state where u 
is parallel to one of the principal axes. This assumption is important because in Section ||| we 
present strong observational evidence that isolated pulsars evolve rapidly towards a stable state of 
this sort, probably under the action of internal friction. Consequently, the properties investigated 
below characterize neutron stars early in their lives. 

The results in this section are mainly presented in the context of a biaxial star (e' = 0) for 
two reasons. First, we find numerically that the rotation properties of biaxial and triaxial stars 
are qualitatively alike, except for pseudo-glitches (Section [T^) which are an intrinsically triaxial 
phenomenon. Second, the aim of this section is to illustrate those aspects of the dynamics that are 
observationally relevant; an exhaustive quantitative survey of the rotation of a triaxial magnet lies 
outside the scope of this paper. Formally speaking, however, biaxiality is a good approximation 
only as long as one has e cr <C e mag (or else e cr S> e mag ) and the internal magnetic field (or crust) is 
symmetric about a unique axis. 



3.1. Precession Period 



The period As of the anharmonic precession is typically ~ |/?| 1//2 |e| 1 , in units of uj 1 . For 
the special case % = 90°, (3 > 0, As is given exactly by 



As 



4 
k 



P 



1/2 



K 



1/2* 



(24) 



where K(k), a complete elliptic integral of the first kind, increases logarithmically from K(0) = ir/2 
to K{1) = 00. Note that the precession is not exactly periodic, because As oc (ij — 7) -1 ^ 2 increases 
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adiabatically on the braking time-scale tq. Indeed, As increases significantly and approaches 
loqTq (so that the separation into slow and fast time-scales in Section |2.4| breaks down) under two 
special sets of circumstances: (i) at k = 1, where u\, u 2 and U3 suddenly swap oscillation modes 
('mode jumping'; see Section [O]); and (ii) at the rj = 7 resonance, where the the star rotates 
steadily with u 2 = U3 = and u% = const for (3 > — 1 (Section |3,4| ). If the star is triaxial, a second 
precession time-scale, As ~ |/3| 1//2 |e'| _1 , is introduced. 

3.2. Precession and Nutation Amplitudes 

Figure |l| displays the precession angle 9, defined to be the angle between u and e3 
(cos# = Us/u), as a function of time. We see from Figure |l] that the star precesses (9 7^ on 
average) and nutates (9 oscillates in a range 9\ < 9 < 9 2 during one precession period), and that 
the slow evolution of the precession and nutation amplitudes is determined by 0, x an d the initial 
orientation of u . In the regime /? > 1 where the free precession period t 2 is shorter than both 
torque-related time-scales To and t\, the nutation amplitude is small and 9 decreases exponentially 
to zero for \ = 20° (thick band in Figure ||). This result, and a similar calculation for \ = 70° (not 
shown), confirm Goldreich's (1970) conclusion that the slow evolution of 9 in the regime (3 ^> 1 
depends solely on x-, with 9 — > 0° for \ < Xcr = cos _1 (3 -1 / 2 ) « 55° and 9 — > 90° for x > Xcr- 

Goldreich's (1970) conclusion is invalid when t\ and t 2 are comparable. The solid and dotted 
curves in Figure |l| both correspond to x = 20° < Xcn but with (3 = 0.8 and hence t% ~ t 2 . Neither 
curve behaves as predicted by Goldreich (1970): either 9 approaches a constant non-zero value 
while the nutation amplitude decreases to zero (solid curve), or else 9 remains constant on average 
with a peak-to-peak nutation amplitude of w 50° (dotted curve). In each case, the precession is 
persistent, and its character is determined by the initial orientation of u}. 

3.3. Do the Magnetic and Rotation Axes Align? 

The magnetic inclination angle a between u and m is defined in terms of ui, u 2 and K3 
by ncosa = u\ sinx + 113 cosx- In Figure |2j, we plot a as a function of time for several choices 
of (3 and x- For x = 0°, (3 = 0.5, we see that a decreases exponentially to zero on the braking 
time-scale. By solving (jl|) to @ analytically for x = 0°, one can show that a approaches zero 
for arbitrary 0, implying that a star subject to a predominantly magnetic deformation ultimately 
becomes an aligned rotator — a state in which it cannot be detected as a pulsar. An aligned final 
state (a = 180°) can also be attained by a star with x 7^ 0° (upper dotted curve in Figure |2|). 
In contrast, if the star nutates persistently as discussed in Section |3.2| , a mimics 9 and oscillates 
within a range (typically tens of degrees). 
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3.4. Oblate and Prolate Stars, and the 77-7 Phase Plane 

An instructive way to view the evolution of the rotation is to follow the star's trajectory 
[i](s), 7(5)] on the 77-7 phase plane. As shown in Section ^4|, 77 and 7 are approximately constant 
over one precession period, varying slowly on the braking time-scale To- Figure |3| shows phase 
diagrams for three stars with x = 90° and different ellipticities. Since 77 and 7 are not exactly 
constant over one precession period, the trajectories are slightly irregular. 

Figure ||(a) shows that an oblate star {(3 > 0) evolves asymptotically towards a state with 
7 = 77 7^ and hence U2 = ^3 = 0, u\ = const (sec and [ 20 1 ) ; in other words, for x = 90°, an 



oblate star always approaches steady-state rotation with u parallel to m (= ei). Equations (|l|) 
to (||) imply that this is a singular fixed point which exists for x = 90° only. The evolution of a 
prolate star (/? < 0) depends on the relative magnitudes of the precession periods t\ (radiative) 
and T2 (free). When the radiative precession is faster (— 1 < j3 < 0), as in Figure |3|(b), the star 
evolves to a state with 7 = 77 7^ 0, as discussed above. When the free precession is faster ((3 < — 1), 
as in Figure ||(c), the star evolves to a state with 7 = 77 = 0, and u does not necessarily align 
with any preferred axis on the way, although it may do so for specific initial conditions. In all 
the above cases, the phase-plane trajectories are confined within the triangular region < 77 < 1, 
min(— 77/3, 0) < 7 < max(77, — 77/3). 

How are the phase diagrams modified for x 90°? Firstly, the trajectories are confined to a 
smaller (larger) triangular region for f3 > (/? < 0) defined by 

< 77 < 1, min(0,7i,7 3 ) < 7 < max(7i,72,7 3 ), (25) 



with 71 = 7] (cos 2 x — 0)i 72 = V s™ 2 X an d 



sin 2^ 

73 = 77[sin 2 (x + V0 -/3sin 2 ^], tan2^ = ^ — . (26) 

(3 - cos 2x 

Secondly, although the (3 > —1 trajectories approach the diagonal line 7 = max(7i, 72, 73), they 
do not stop there like in Figure ||. Instead, they bend downwards to merge with the diagonal, 
travelling down along it as 77 decreases. This is because the state 7 = max(7i, 72, 73) is only a 
fixed point for x = 90°, as pointed out above. 

The trajectories in Figures |3|(a)-(c) asymptotically approach, or travel exactly along, the 
lines 7 = 77 and 7 = —77/? except in the case — 1 < (3 < (Figure ||[b]), where trajectories with 
7<— (3 at s = subsequently cross the line 7 = —r]/3. When this happens, a phenomenon we call 
'mode jumping' takes place. As shown in Figure ||, u\ and U2 interchange oscillation modes at 
s ~ 2.7 u>oTq, swapping between a wine-glass mode and a sinusoidal mode, while the 113 oscillation 
becomes temporarily flatter-peaked. 



-12- 



3.5. Pseudo-Glitches in the Frequency Derivative of a Triaxial Star 

A triaxial ellipsoid of inertia arises naturally if, for example, the deformation is predominantly 
magnetic with appreciable quadrupolar and off-centred components, or if one has e cr ~ e mag . 
Figure || plots the angular frequency derivative u as a function of time for a triaxial star with 
e' = 0.09e. We see that the smooth, braking-related decrease of \u\ is punctuated by sudden, 
quasi-periodic spikes in which \u\ changes by up to 90 per cent. We call these excursions 
'pseudo-glitches'. They resemble true glitches because (i) they recur quasi-periodically with period 
~ xqTq 10 yr for the Crab), and (ii) u returns to its trend value after each excursion. However, 
they are manifestly not true glitches because (i) their rise time is too long (cf. <J 10~ 10 to for Crab 
glitches), and (ii) they do not cause u itself to increase. Note that ii oscillates about the average 
spin-down trend even for e' = 0, but not in the spiky fashion of Figure ||. 

What is the physical origin of pseudo-glitches? The dotted curve in Figure [5] shows that 
pseudo-glitches coincide with rapid changes in a, accompanied by mode jumping, which occur 
when the free precession is modulated on a fast time-scale (e'/e)r2 which couples resonantly to the 
radiative precession (i.e. eri i=a e'r2). It turns out that the sharpness of the spikes is sensitive to e' 
and uG(xo); for the example in Figure ||, the spikes are washed out once e' falls outside the range 
0.05 <; e'/e 0.15. Figure |5| shows the case x = 40°. For larger x, the small bumps between the 
spikes increase in amplitude until they become spiky themselves. For smaller x, the bumps flatten 
until they disappear. 

3.6. Corotating Magnetosphere 

The foregoing results pertain to the regime xo <S 1, where ro is taken to be the stellar 
radius R. However, recent work by Melatos (1997) suggests that the corotating magnetosphere 
of a neutron star acts as a perfectly conducting, rigid extension of the stellar interior out to 
a characteristic radius r v where outflowing plasma is not constrained to flow along magnetic 
field lines by cyclotron losses, and that it is therefore necessary to set ro = r v when calculating 
electromagnetic spin-down properties like pulsar braking indices. For young pulsars (e.g. the Crab, 
PSR B1509-58, PSR B0540-69), one finds r v £ c/uj, and hence x £ 1. 

In the regime xq <^ 1, the radiative precession period satisfies t\ ~ to, the near-field radiation 
torque uncouples from the free precession, and the evolution of the precession amplitude is 
governed completely by x/Xcr (Goldreich 1970). However, a pulsar born with xq ^ 1 soon evolves 
towards the regime xo <C 1 as to decreases. Therefore, unless the free precession amplitude 
decreases to zero during the initial braking phase with xq <J 1 {not the outcome in general), the 
rotation behaves thereafter in the way described in Sections [O] to [Q| . 
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4. APPLICATION TO OBSERVATIONS 

Timing and polarization studies suggest that most, if not all, young pulsars do not precess in 
the manner described in Sections ^ and ||. The implication is that isolated neutron stars born with 
a large precessional motion approach stable spin states (u) parallel to a principal axis of inertia) 
over times that are short compared to their current ages, probably due to internal friction. We 
outline the conditions for a neutron star to be born with a large precessional motion in Section 



4.1, summarize the status of observational searches for precession in Section |4.2|, and discuss the 



observational consequences of frictional stabilization in Section 4.3 



4.1. Precession at Birth 

When a neutron star is born, it spins about an axis u»o dictated by conservation of angular 
momentum (of the degenerate remnant and the ejecta) during the supernova explosion. There is 
no reason why wo should immediately be parallel to the principal axis of greatest non- hydrostatic 
moment of inertia (the magnetic axis m, since the star is a fluid). Indeed, Thompson & Duncan 
(1993) argued that post-collapse convection destroys any correlation between ujq and m. 

Viscous dissipation in the fluid star forces ujq to approach m over time. The dissipation rate 
is therefore critical in determining whether the star will exhibit a large initial precession when its 
crust crystallizes. Unfortunately, the viscosity of a newly born neutron star is poorly constrained. 
Cutler & Lindblom (1987) estimate the viscous damping time (e.g. for stellar oscillations) to be 
roughly 3 x 10 2 (p/10 17 kgm~ 3 )~ 125 (T/10 9 K) 2 yr, to be compared with a crust crystallization time 
of less than one year, but this estimate is known to be valid only in narrow ranges of density p 
and temperature T centred on the above fiducial values. In what follows, we make no assumption 
about the viscous damping time and instead explore several possible scenarios. 

If the viscosity is high enough, wo aligns with m first, before the crust crystallizes. Assuming 
that the symmetry axis of the crust when it crystallizes is along ujq (likely, though not certain), 
then the principal axis e3 (from both elastic and magnetic contributions) is parallel to uq, and 
there is no precession. 

If the viscosity is low enough, the crust crystallizes first, before wo has time to align with m 

- the order of events implied by the viscosity estimates of Cutler & Lindblom (1987). In this 
scenario, two things can happen: (i) one has e cr 3> e mag , and the symmetry axis of the crust when 
it crystallizes is along uq, so that e3 is parallel to wo and there is no precession; or (ii) one has 
£ cr ^ £ mag^ go fa&t e% is parallel to m — which is not parallel to wo at the epoch of crystallization 

— and there is a large precession. Below we explore the observational implications of scenario (ii). 
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4.2. Changes in Pulse Profile and Polarization 

The results of Section ^ imply that a neutron star in an arbitrary initial spin state precesses 
(and nutates) with period ~ xqtq and a typical amplitude of tens of degrees, and that the motion 
is persistent in general. One therefore expects fractional changes of order unity in the pulse profile 
(e.g. relative height or separation of conal components), magnetic inclination angle a (measured 
from polarization-swing data), and angular frequency derivative u (measured by timing) over a 
single precession period. For young, Crab-like objects (xo 10~ 2 , tq 10 3 yr), the fractional 
changes amount to ~ 10 per cent per year and ought to be readily observable; for old objects 
(xo ~ 10~ 3 , To ~ 10 6 yr), the changes amount to ~ 0.1 per cent per year and are harder to detect. 

Contrary to the above expectation, the only reliable detection of neutron-star precession to 
date has been Weisberg et al.'s (1989) discovery of the general relativistic geodetic precession 
of PSR 1913+16. Six years of accurate measurements of the doubly-peaked radio pulse profile 
revealed that the flux density in the first component is decreasing relative to the second by ~ 1 
per cent per year, consistent with the line of sight moving across a spot of enhanced emissivity 
in the magnetosphere at the rate prescribed by geodetic precession. Weisberg et al.'s (1989) 
upper limit [Aio| < 0°.06 on the six- year change in pulse width w implies a maximum precession 
amplitude <J 0°.4yr , consistent with a relativistic origin and much larger than the predicted 
radiative precession amplitude of ps 0°.001yr _1 . Hence PSR B1913+16 does not usefully constrain 
the amplitude of radiative precession in old pulsars. 

No unambiguous instances of non-relativistic precession are known. Lyne et al. (1988) 
reported a quasi-sinusoidal variation with a period of ~ 20 months in the phase residuals of six 
years of Crab timing data, and a similar feature was claimed to exist in Vela (McCulloch et al. 
1990). However, the variation may be an artifact of an unexpectedly high Co during the exponential 
relaxation following an overlooked glitch. The quasi-periodic nature of the Vela glitches (Lyne et 
al. 1996) is also suggestive of precession, but the observed wasa function of time does not resemble 
Figure |5[ Ulmer (1994) reported that the intensity ratio of the two peaks in the gamma-ray pulse 
of the Crab seems to vary sinusoidally with a period of ~ 14 yr in both the 50-500 keV and 50 MeV 
bands (but not in the optical). However, the 14-yr period differs from that observed by Lyne et 
al. (1988) and a 60s modulation of the Crab's optical pulses found by Cadez et al. (1997). 

Direct measurements of a from pulsar polarization swings do not show any evidence for 
secular changes on a yearly time-scale. It has been claimed that changes in a on a time-scale 
of ~ 10 4 yr explain the braking indices of the Crab and Vela, and sudden jumps in a explain 
the persistent increase in uj following a glitch (Allen &: Horvath 1997; Link & Epstein 1997), 
but the relevant time-scales are too long and too short, respectively, to be precession-related. 
Tauris & Manchester (1997) used polarization-swing data for more than 100 pulsars to construct 
the a distribution of the pulsar population, corrected for decreasing beam radius with age. The 
corrected distribution is skewed towards small a, with (a) decreasing on a time-scale of ~ 10 7 yr 
— once again, too gradual to be a precession-related effect. 
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The observational evidence against substantial changes in pulse profiles and polarization 
properties implies that any precession initially present is damped rapidly. Consequently, direct 
detection of a precessing, isolated pulsar may only be possible in the immediate aftermath of a 
Galactic supernova. If the precession amplitude is sufficiently large, the newly born pulsar will 
shine intermittently, as its emission cone drifts into and out of the line of sight, as well as exhibit 
the profile and polarization changes discussed above. 

4.3. Crustal Heating by Internal Friction 

Frictional damping inside a neutron star proceeds rapidly in theory. Several authors have 
estimated the dissipation times due to time-dependent elastic strain (Chau & Henriksen 1971; 
Macy 1974) and imperfect core-crust coupling (Shaham 1977; de Campli 1980; Alpar & Ogelman 
1987; Link et al. 1993) and found them to be small compared to the present ages of young, 
Crab-like pulsars. (The estimates assume a small precession amplitude; cf. Section ||.) In addition, 
there is the analogy of friction inside the Earth. Dynamic satellite measurements have revealed 
that the Earth's non-hydrostatic ellipsoid of inertia has its e3 axis parallel to u to an excellent 
approximation^ (Lambeck 1980, p. 31), and that the direction of e3 fluctuates by less than l" 
under the action of solar and lunar tides, compared to ~ 10° for the ei and e 2 axes (Bursa & Pec 
1993, p. 227). Clearly, any drift of u away from e3 is rapidly damped. 

In this paper, we are not concerned with the precise origin of the friction in a neutron star; 
we simply suppose it exists and examine the fate of the dissipated energy. Initially, when the star 
is precessing, its angular momentum and energy are given by L\ = [if (u^ + uj^i) + /fwfj] 1 / 2 and 
E{ = l/i(o;f i + u>2i) + ilswfi. After the precession has been damped, the final angular momentum 
and energy are given by Lf = I^ui^ and Ef = ^i^ff (u — ► e3 by analogy with the Earth). If the 
damping occurs over a time that is short compared to tq and t%, we have L; ~ Lf and hence 

AE = -\eh{ujl+uj 2 2 {). (27) 

For wii ~ ~ u>o (Section ||), the total dissipated energy is of order eI±UQ. 

We now suppose that the dissipation is localized in the inner and outer crusts, where 
the ions are organized into a lattice and the shear modulus is non-zero, and we estimate the 
resulting increase in the star's thermal luminosity L. The thermal conduction time in the crust, 
T con d, is not known with certainty (Nomoto & Tsuruta 1987), so we appeal to the extreme 
cases of slow cooling (ra <C T con d) and fast cooling (r^ 3> T con d) to place bounds on L. In 
the regime T con d, the dissipated energy AE is stored in the crust as heat for a time 

Tcond- The crustal heat capacity at a temperature T is given by c v m 312k B (T/e D ) 3 (Shapiro & 
Teukolsky 1983, p. 100) in the regime where the Debye temperature 8 D ~ 2 x 10 10 K satisfies 



5 1 thank P. Goldreich for bringing this fact to my attention. 
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T <C D (i.e. at a density p ~ 1 x 10 14 gcm -3 ). Taking AE 1 = eliLO^, we find that the final 
temperature of an iron crust of mass 10 -3 M is Tf = 2 x 10 9 (ecjQ) 1//4 K, and its thermal luminosity 
is L = AirFPaT^ = 1 x 10 46 eu;o ergs -1 . In the opposite regime S> t coi1c j, the dissipated heat is 
conducted rapidly through the crust and one has L = AE/rd = 1 x W^cujqT^ 1 ergs -1 (cf. de 
Campli 1980, p. 308), less than the former value for Td > 0.1s. 

The above estimates imply that the minimum thermal luminosity of a newly born pulsar due 
to frictional damping of its radiative precession is 

Lf»3x 10 31 ( ( " X ( ^-Y 1 ergs -1 . (28) 
V10 -12 / V10 3 rads -1 / \lyr/ 

This ought to be detectable given Td <^ 10 -3 yr even if there is significant magnetospheric X-ray 
emission beamed towards the observer and with Td <J 3yr if there is not. The duration of the 
thermal source is the maximum of Td and r con d, starting from the time when the star first 
crystallizes into a body with a rigid crust spinning down electromagnetically. If the crystallization 
epoch occurs very shortly after the supernova explosion itself, it may not be possible to detect L 
at all. 



5. CONCLUSION 

In this paper, the rotation of a rigid, aspherical, internally frictionless neutron star is analyzed. 
We show that, in general, the free precession period T2 due to elastic and magnetic deformations 
is comparable to the radiative precession period t\ associated with the near-field component of 
the radiation torque. In other words, the 'free' precession is not truly free, a fact that has been 
overlooked in the literature to date. 

In the regime t\ ~ T2, the star rotates in a distinctive way: (i) it precesses and nutates 
anharmonically, typically with an amplitude of tens of degrees (Section |3.2[ ); (ii) the magnetic 
inclination angle a swings through tens of degrees during one precession period (Section [3l|); (iii) 
the precession can persist or decay to zero (i.e. steady rotation) depending on the parameters 
j3 and x an d the initial orientation of u (Section |3.4| ); and (iv) the frequency derivative u < 
oscillates about its overall spin-down trend, exhibiting spiky, glitch-like behaviour for triaxial stars 
with €T% « e'r2 (Section |3.5|). 



The precession and nutation lead to fractional changes of order unity in the pulse profile, 
polarization swing and u of an isolated pulsar on a time-scale ~ xqTq <C tq, with xqTq ~ 10 yr for 
young, Crab-like objects and xqTq ~ 10 3 yr for old pulsars. Such changes are not observed. One 
plausible explanation is that a young neutron star has u parallel to e3 when its crust crystallizes 
shortly after birth — but this is not true for neutron stars with large hydromagnetic deformations, 



given current viscosity estimates (Section 4.1). Another explanation is that the precession and 
nutation are damped by internal friction, perhaps due to time-dependent elastic strains in the 
crust. If the damping takes place over a time Td, we show (Section |4,3j ) that the dissipated energy 
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AE « cIiujq either heats the crust to a temperature Tf = 2 x lO^eu^) 1 / 4 K for ra <C T con d, 
yielding a thermal X-ray luminosity L = 1 x 10 46 eu;g erg s _1 , or else is conducted rapidly to the 
surface, yielding L = 1 x 10 45 ea;Qr d ~ 1 erg s _1 . The luminosity L may be detectable depending on 
how soon after the supernova explosion the neutron-star crust crystallizes, and the intensity of 
magnetospheric X-ray emission at that epoch. 

I thank Sterl Phinney, Peter Goldreich and Jon Arons for comments. This work was supported 
by NASA Grant NAG5-2756, NSF Grant AST-93-15455, and the Miller Institute for Basic 
Research in Science. 
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A. DERIVATION OF EULER'S EQUATIONS 



We ignore the slight distortion of the star from a perfectly spherical figure when calculating its 
radiation fields and the radiation-reaction torque. The electromagnetic fields E(x, t) and B(x, t) 
generated by a magnetized, conducting sphere rotating in vacuo were derived using a multipole 
method by Deutsch (1955) and subsequently corrected for minor typographical errors by Melatos 
(1997); see also Kaburaki (1981). The radiation torque exerted on the rotating sphere can be 
calculated from E(x, t) and B(x, t) by integrating the angular momentum flux vector over any 
surface S enclosing the sphere: 



N = e / [(x x E) E • dS + c 2 (x x B) B • dS - \{E 2 + c 2 B 2 )(x x dS) 



(Al) 



Let the radius of the sphere be ro, let its angular frequency be u, let a denote the angle between 
its rotation and magnetic axes, and assume that the frozen-in magnetic field is dipolar, with polar 
magnitude Bq. Then, in a Cartesian coordinate system (x,y,z) where the z axis is oriented along 
the instantaneous rotation axis u and the magnetic axis m simultaneously lies in the x-z plane, 
the instantaneous radiation torque assumes the form 



(N x ,N y ,N z 



3 r 



with 



F(x ) 
G(x ) = 



sin a cos aF(xo), sin a cos aG(xo), — sin 2 aF(xo) 



5{x% - 3x% + 36) 
3(x§ + 6) 



5x (x(j 



+ 



+ 



3(xg + l)' 
3 - 2x 2 Q 



(A2) 



(A3) 



(A4) 



■-, -3xg + 36) 15x (x5 + l) 

Strictly speaking, it is incorrect to refer to (A2) as an instantaneous torque, because it is calculated 
under the assumption that the star has been rotating, and will continue to rotate, at a constant 
angular frequency uj; formally, it is assumed that the radiation fields exist for all t and are 
proportional to e tujt , as for an infinitely massive star. However, the approximation is an excellent 
one for a neutron star, where the braking and precession time-scales To, t%, T2 (Sections |2~ 
satisfy r , n, r 2 » a; -1 . 



The form factors (A3) and (|A4j ) differ from the expressions F(xq) = 1/3 and G{xq) = 1/2xq 
appearing in previous works (Davis & Goldstein 1970; Goldreich 1970). There are two reasons for 
this difference: (i) previous authors only included terms of leading order in the small parameter 
xo (cf. Section |3.6| ), whereas ( |A3| ) and (A4) are exact for arbitrary xq; and (ii) previous authors 
modelled the star as a sphere of uniform internal magnetization, whereas in this paper the star 
is modelled as a perfectly conducting sphere with a point magnetic dipole at its centre (Deutsch 
1955), affecting the polynomial coefficients in ( |A3[ ) and (|A4j ). 

The Cartesian coordinate system (x,y,z) rotates with respect to the star as cj changes 
orientation under the action of N. Among other things, this causes a to change with time. In 
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order to write down Euler's equations, it is necessary to reexpress N in body coordinates that 
are fixed with respect to the star. We choose the body axes to be the principal axes of the star's 
ellipsoid of inertia, ei, e 2 and 63; the magnetic axis m, which is also fixed with respect to the star, 
is taken to lie in the plane spanned by ei and e3, at an angle x to e3 (a slight loss of generality). 
The transformation from (x, y, z) to body coordinates is time-dependent, but it is not an Euler 
transformation because (x,y,z) is a non-inertial frame. It is described by a matrix [Aij], defined 
through (ei,e2,e3) = [Aij] • (i,j,k), which can be represented in the form 

cos <f> cos 9 cos ip — sin <p sin ip sin <p cos 9 cos ip + cos <p> sin ip — sin 9 cos ip 
[Aij] = j — cos <p cos 9 sin ip — sin <p cos ip — sin (p cos 9 sin ip + cos (p cos ip sin 9 sin ip | , (A5) 
cos (p sin 9 sin <p sin 9 cos 9 

where the angles (p, ip, 9 (not Euler angles) depend on t through the principal angular velocity 
components U\, uj 2 , W3 as follows: 

lo cos x- ^3 cos a 

cos(p = — j—, (A6) 

sm a{uj z — uj^) l ' z 

. u) 2 sinx 

smcf) = — : — jus, (A7) 

sm a(uj z — > 

cosip = — Ul 2n/2 , (A8) 

Sil1 ^ = T~2 2U72"' ( A9 ) 

cos 9 = (A10) 
sin0 = [1-ujI/uj 2 ) 1 / 2 . (All) 

The angle a is chosen to lie in the range < a < ir and satisfies 

lo cos a = oj\ sin x + <^3 cos x- (A12) 



Upon substituting ( |A6| ) to ( A12j ) into (A5) and employing the relation (Ni, N 2 , N3 



[A^] ■ (N x , N y , N z ), we arrive at the principal components of the radiation torque featured on the 
right-hand sides of Euler's equations ([[]), (||) and 



B. APPROXIMATE ANALYTIC SOLUTION FOR e ' = 0, x = 90° 

In this appendix, we derive an approximate solution to Euler's equations ([[]), (Q) and (||) for 
a biaxial star (e' = 0) with x = 90°. The solution is accurate provided the precession periods n 
(radiative) and t 2 (free) are much smaller than the braking time-scale To- 

When the slow braking terms, proportional to (woTo) -1 F (xo) , are neglected relative to the 
fast precessive terms, proportional to e and (ujqTq)~ 1 uG(xq), Euler's equations reduce to the 



-21 - 



zeroth-order system (|i~5| ) to ([P?]) which, in the special case \ = 90°, reduces to 

ui = -eu 2 u 3 , (Bl) 
ii 2 = e{l + (3- 1 )u l u z , (B2) 
U3 = -e0~ 1 uiu 2 . (B3) 

The solutions of ( |B1| ) to (B3) are Jacobian elliptic functions. The relative signs of the coefficients 
on the right-hand sides determine the solution branch, given the physical requirement that u±, 
U2 and 1*3 are real quantities whose squared amplitudes are non-negative. We distinguish three 
solution branches: 



Case I: > 0. 



Ml 

U 3 
k 2 
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i + v0\ 1/2 
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sn 



Case II: -1< < 0. 
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Case III: < -1. 
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In (B4) to ( |B15| ), the quantities 77, 7 and 3> are all constants of the motion (see Section 2.4). 



The trivial cases (3 = — 1 (harmonic precession about and (3 = (spherical star; see Davis &: 
Goldstein 1970) are not treated here. 

When terms proportional to (ujqTq)~ 1 F (xq) are restored to Euler's equations, the above 
solutions remain approximately valid, but the constants of the motion are converted into slowly 
varying functions of s. We now compute the slow evolution of r/(s), 7(5) and ^(s). First, we 
substitute ([B4| ) to ( B15|) into ([[]) , © and (||) and perform the time derivatives explicitly. For each 



solution branch, this results in a system of three equations linear in 77, 7 and <3?. For example, the 
[3 > branch yields 

7 + D/3 _ F(xo) Vm9t (B16) 



2(7 + Tj/3) UJqTq 

13 \ 1/2 r)-7 



-s + $ 



dn2 = _F(xo>? sncndn) (B17) 



j + 2 , i)-7 ,2 ^(^o)^ , 2 , e 2 (7 + , 
-T7 -cn H ; dn = —- r dn +- ; ^ ln , -— 7-sncndn. (Bio 

2/3(77-7) 2(77-7) (l + e)u T {i + e )^/2^_ 1 y/2 V ) 

In ( JB16| ) to ( piS| ), Euler's equations are linearly combined in such a way that the coefficients of 
7), 7 and $ are squares of elliptic functions, in order to ensure a non-trivial result when averaging 
over The (omitted) arguments of the elliptic functions are the same as in (|B4j ) to ( |B6[) . We 
now average over holding constant 77, 7, and their derivatives. The results for the three solution 
branches are as follows: 



Case I: (3 > 0. 



2F(xoM(l + e)( 7 + Vf3)h(h + /H3) + - 7K3] (mg) 



u> T {l + e)P{(3 + l)h 
2F(xq)t/[(1 + e)( 7 + r?/?)Ji(/ 2 - J 3 ) + P(r] - l)l 
u T (l + e)((3 + l)I 3 



(B20) 



* = ^WTWh W^J es - (B21) 



Case II: -1< (3 < 0. 

2F(x )t7{(1 + e)/3(77 - 7 )J 3 + 7 J 1 [/3J 3 - {£ + Ijjg]} 
a; r (l + e)/3I 3 



(B22) 



2F(x )r n I 1 
u; t (1 + e) 

i ^_^/r/- 7 V /2 F(x )77[(l + e)/?(77-7)/ 3 -(/? + l)7/i/ 2 ] / -fl V /2 

* = (TT^3 1— J ^0(1 + ^/3 I^FtV 65 ' (B24) 



Case III: (3 < -I. 



2F(x ) V {(l + e)( 7 + V f3)I 2 [(f3 + 1)J 3 + 7j] - Q3 + 1)7/3} 
WO To(l + e)/3(/? + l)/ 3 



(B25) 
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2F(x )r l [-(l + e)( 7 + vP)hh + (/? + 1)7/3] (B2g) 



w ro(l + e)(/3 + l)/ ; 



* = Mi + ^T^ (B27) 



In ( |B19| ) to (B27), the <E> averages I{ (1 < i < 4) are all functions of the modulus A; appropriate for 
each branch and are defined by Ii = (sn 2 <3?), I2 = (cn 2 <£), I3 = (dn 2 <I>), and Z4 = (cn 2 <£ dn 2 <£). 

In the limit j3 — ► 00 (i.e. A; — ► 0), the three solution branches merge together into one, the 
elliptic functions reduce to trigonometric functions, and the equations for 77, 7, $ reduce to those 
given by Goldreich (1970). 
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Fig. 1. — Precession angle 6, in degrees, as a function of time, in units of the spin-down time-scale 
r . Thick band: (3 = 20, x = 20°, uG(x ) = 47.4, u 1>0 = 0.387, u 2 ,o = 0.6, u 3fi = 0.7. Although 
not apparent to the eye, the band is a rapid oscillation with a peak-to-peak amplitude of 1.6°. 
Solid curve: j3 = 0.8, x = 20°, uG(x ) = 47.4, u lfi = 0.387, u 2 , = 0.6, n 3i0 = 0.7. Dotted curve: 
[3 = 0.8, X = 20°, uG(x ) = 47.4, u 1>0 = 0.387, u 2 ,o = 0.6, u 3 ,o = -0.7. 
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Fig. 2. — Magnetic inclination angle a, in degrees, as a function of time, in units of the spin-down 
time-scale tq, for uG(xq) = 47.4 and initial conditions u± t o = —0.5, U2 t o = 0.4, U3 5 o = 0.768. Lower 
solid curve: (5 = 0.5, x = 0°. Lower dotted curve: (5 = 0.5, x = 45°. Upper solid curve: /3 = —0.5, 
x = 60°. Upper dotted curve: (3 = -0.5, x = 90°. 
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Fig. 3. — Phase diagrams showing trajectories in the 77-7 plane, for uG(xq) = 47.4, x = 90°, and 
< s < 10 2 woTo. The trajectories begin at 77 = 1, with different initial values of 7 corresponding to 
different initial conditions, and move from right to left across the diagrams. The broken diagonals 
are the lines 77 = 7 and rj = Phase diagrams are shown for three stars: (a) (3 = 0.5 (oblate), 

(b) (3 = -0.5 (prolate), (c) (3 = -1.5 (prolate). 
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Fig. 4. — Principal angular frequency components u\ (dashed curve), 112 (dotted curve) and 113 
(solid curve), in units of loq, as functions of time, in units of the spin-down time-scale To- Note 
that u\ and U2 interchange oscillation modes at s « 2.7luqTq. All curves are for \ = 90°, (3 = —0.5, 
uG(x ) = 47.4, u 1>0 = 0.539, u 2fi = 0.7, u 3j0 = 0.469. 
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Fig. 5. — Angular frequency derivative ii (solid curve), in units of {ujqtq)~ 1 , and magnetic 
inclination angle a (dotted curve), in degrees, as functions of time, in units of the braking time-scale 
To- The local minima and maxima of u and a correspond one-to-one. Both curves are for x = 40°, 
P = 9.5, e ' = 0.09e, uG(x ) = 9.8, u lfi = 0.873, u 2 ,o = 0.28, u 3fi = 0.4. 



